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1. Introduction 


1 . 1 General 

The safety and reliability of structures has always been a 
matter of vital concern to the aerospace industry. In this 
respect, fracture mechanics (FM) is a specially useful 
technology, since it can provide a quantitative description of 
the capability of structural parts to tolerate flaws. Initially, 
FM concepts covered quasi-linear elastic conditions (LEFM) . 

Later, these methods were further developed to cover more general 
situations. Specifically, there was a need to extend these 
concepts to include cases where yielding was not necessarily 
contained in very small regions, for the case of new and tougher 
materials, higher loads, thinner sections, et cetera. This led 
to the development of the so-called Elastic Plastic Fracture 
Mechanics (EPFM) Methodology. 

To apply these methods, two pieces of information are 
needed: the so-called material/specimen response to deformation, 

and the material response to crack extension. The former, 
obtained by finite element analysis or experimental calibration, 
consists of two expressions connecting the J-integral, load P, 
load-point displacement v, and crack length a for the specimen 
geometry of interest; the latter consists of a characterization 
of the way the material resists crack extension for the type of 
load applied: J, or a similar parameter, versus crack extension, 

for monotonic load, da/dN versus AK or AJ for cyclic loading, 
da/dt versus K, C* or C t for creep crack growth, et cetera. It 
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is assumed that within some limitations, these curves are 
specimen geometry independent, i.e., the curve obtained from a 
small laboratory specimen applies to the structural part under 
consideration, as well. 

A simple computer program can be developed to combine the 
two pieces of information mentioned and assess the structural 
reliability of the structural part of interest. 

It is very important to devote effort to guarantee that the 
curve of material response to crack extension is, in fact, 
geometry independent. That is, it is important to understand the 
limitations of the parameters and/or approaches used, identify 
clearly their limits of validity, and eventually improve the 
characterization of the phenomenon, proposing new parameters and 
methods to extend the range of applicability of existing models. 

1.2 Elastic Plastic Fracture Mechanics 

Specifically, for the case of EPFM applied to monotonic 
load, the mentioned limitations are expressed in terms of the 
amount of crack extension to ligament ratio, r, the ratio of 
ligament to applied J over the yield strength, m, and the ratio 
of logarithmic increase of J to logarithmic decrease in ligament, 
0 ). 

To overcome some of these limitations, particularly the one 
on r, Ernst [1.1] proposed a modified version of J called J M . 
Resistance (R) curves plotted in terms of J M were not subjected 
to the same limitations as those using J, and in general, showed 
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a better correlation between specimens of different size and 
geometry. 

More recently, this methodology was further extended: 
general formulas were developed for and Jp for growing cracks, 
criteria were proposed to identify the limits of applicability of 
both parameters, methods were presented to make use of the 
information of experimental points beyond this limit, and several 
schemes were proposed to extrapolate small laboratory specimen 
resistance curves to large amounts of crack extension, using J M , 
Jq, or other parameters [1.2]. 

Although the progress made has been significant, and 
understanding has been gained on how to represent the R curve 
[1.2-1. 6], there are still several very important points that 
need to be addressed before the method can be safely applied. 

Among the most important ones, the need to extend this whole 
methodology to include cases involving three dimensions (3D) must 
be mentioned. 

Specifically, it is mandatory to know how specimen 
thickness, constraint, and the possible dependence of the 
fracture mechanism on specimen thickness may affect the fracture 
resistance . 

Ultimately, this knowledge gained from "2D" planar specimens 
should be used to explain and predict the behavior of real life 
3D defects found in structures, i.e., surface or embedded cracks, 


et cetera. 
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1.3 The Leak Before Burst (LBB) Criterion 

Pressure vessels containing surface flaws are often required 
to comply with the so-called LBB criterion. LBB is understood as 
the condition in which an assumed initial flaw will grow through 
the wall of a pressure vessel and cause leakage rather than 
bursting . 

In particular, pressure vessels of interest to NASA have to 
comply to MIL-STD-1522A Standard General Requirements for Safe 
Design and Operation of Pressurized Missile and Space Systems. 
This document requires that: (1) a/2c (crack depth to total 

width ratio) needs to be in a range from 0.05 and 0.5, and (2) 

LBB will occur if Kic/°op > 2aB ( - ) *^ with aa op < Oy S and a > 1. 

The rationale behind this expression is that the initial 
semi-elliptical flaw will grow in a self-similar manner, i.e., 
keeping a/2c constant until the crack depth a is exactly equal to 
the thickness B, as shown in Figure 1.1. At that time, it is 
considered that the flaw becomes a through crack with a total 
length of 2c, with 2c = (2c/a) 0 B, where the subscript 'o' stands 
for 'initial'. Finally, to prevent the crack from running 
unstably in the longitudinal direction, i.e., bursting, it is 
then required that the toughness Kj c be bigger than the applied K 
given by K app = a op (rcc) 0 - 5 . 

The weaknesses of this Standard are these: (1) the above 

equation only holds for a/2c = 0.5, (2) the flaw shape is 

considered to always remain elliptical with constant a/2 c, and 
(3) the whole analysis is based on LEFM concepts. 
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On the other hand, for the real materials, thicknesses, and 
typical flaws of interest, the situation is markedly different 
from the one assumed above, as can be seen in Figure 1.2. The 
cracks, clearly, do not grow in an elliptical self-similar 
manner, but rather in a very complex shape, with a dimension in 
the direction parallel to the surface, longer in the interior 
than on the surface. Moreover there is no guarantee that this 
dimension can be conservatively estimated by taking the 
original (2c/a) o and multiplying by B. 

1.4 This Project 

The EPFM Methodology has evolved significantly in the last 
several years. Nevertheless, some of these concepts need to be 
extended further before the whole methodology can be safely 
applied to structural parts. Specifically, there is a need to 
include the effect of constraint in the characterization of 
material resistance to crack growth and also to extend these 
methods to the case of 3D defects. 

As a consequence, this project was started as a 36 month 
research program with the general objective of developing an 
elastic plastic fracture mechanics methodology to assess the 
structural reliability of pressure vessels and other parts of 
interest to NASA containing defects. 

The project is divided into the following tasks. 
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Task 1. Constraint and Thickness Effects 

This task includes the study of the problem of 
constraint and thickness effects, in different specimen 
sizes and geometries in materials of interest. 

Specifically, the following subtasks will be performed: 

a) The large body of available data from centers 
around the World will be gathered to study this 
effect in specimens of different size and 
geometry . 

b) Resistance to crack growth tests will be conducted 
using specimens of different size and geometry, on 
at least one material of interest. The material 
will be provided by NASA; Georgia Tech will 
machine the lab specimens. 

c) Characterization of fracture surfaces to determine 
mechanisms of fracture, and typical surface 
dimensions will be performed using modern 
quantitative metallographic techniques. 

d) Using the information obtained, models will be 
developed to describe the effect of constraint on 
the growth of cracks under elastic plastic 
conditions . 
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Task 2 . T2ix.ee Dimensional Cracks 

The problem of applicability of EPFM concepts to 3D crack 
problems, in materials of interest, will be studied in this task. 
Specifically, the following subtasks will be performed: 

a) Plates containing surface cracks with different 
initial crack aspect ratios and relative crack to 
plate geometry dimension will be tested. The 
evolution of the crack shape (planar) and the 
crack surface displacement with loading will be 
determined . 

b) Analytical and numerical efforts will be devoted 
to determine values of J and constraint along the 
crack front. 

c) The models and information obtained from Task 1 
will be used here to predict the behavior of these 
3D cracks . 

d) Predictions and experimental resullts will be 
compared and, if necessary, refinement of the 
models will be made. 
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Task 3 . Leak Before (LPB 1 Criteri on 

Finally, the body of information obtained in the previous 
tasks will be organized in a Methodology format to assess the 
structural integrity of parts containing defects, in the spirit 
of the current LBB criterion. 

1.5 This Report 

This report covers the activities of the period March 1993 
through August 1993. In this period, full advantage was taken 
from the experience and knowledge gained in previous projects 
[1.6-1. 9]. In particular, some efforts were devoted in this 
project to complete and extend previously obtained results. 

The report is organized as follows: In Chapter 2, a 

computer modelling algorithm used to simulate the growth of a 
semi-elliptical surface crack is explained in detail. This is an 
excerpt from the thesis of D.W. Boatwright [1.9]. 

In Chapter 3, a finite element investigation is presented. 
This investigation, an excerpt of the thesis of W.J. Curtin 
[1.10], compared the theoretical (HRR) stress field to that 
produced by elastic and elastic-plastic models. The difference 
in these stress fields is the constraint effect. 

In Chapter 4, experimental efforts to characterize three 
dimensional aspects of fracture present in "two dimensional", or 
planar configuration specimens have been continued. This 
discussion specifically contains a preliminary discussion 
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associated with the determination of, and use of, crack face 
separation data. 
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Figure 1 Simple Leak Criterion 
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Figure 2 Growth of a Part-Through Crack to Critical Size 



CHAPTER II: An Excerpt from DW Boatwright's Thesis 


CHAPTER III 


COMPUTER MODELLING OF SURFACE CRACK GROWTH 


Background 

The best way to determine the physical response of a 
cracked body to a particular set of loading conditions is to 
run a test on the structure in question. Unfortunately, this 
is not a practical solution in many cases. One alternative to 
testing is to use a mathematical model that can predict the 
behavior of a test specimen in any given situation. 

In order to use a mathematical model, it is first 
necessary to run a series of mechanical tests to characterize 
to behavior of the specimens. Once these tests have been run, 
it is then possible to use the test data to develop a valid 
model for the specimen in question. Today, libraries of 
functions exist that can be used to model most common fracture 
specimens given the material properties [ 18 , 22 ]. 

After developing a mathematical model that accurately 
predicts test results, the key parameters of the test can be 
changed to observe their influence on test results. 



Currently, the only way to model a surface crack specimen 
is to run a three-dimensional finite element model. This 
modelling process is computationally intensive and very time 
consuming since a new finite element mesh must be constructed 
for each iteration in the solution process. 

It would be far more efficient to develop a special 
purpose fracture mechanics based computer program designed to 
model surface crack behavior. One of the goals of this work 
was to do just that. This chapter covers the development and 
testing of a fracture mechanics based computer program for 
modelling the behavior of surface cracks. 

The intent in the development of this computer program 
was to develop a program general enough to model any surface 
crack geometry in any material that is under J-dominant 
conditions. Given the material deformation and fracture 
properties, the program should accurately predict the 
distribution of crack growth along the front. 

The program described in this chapter represents a 
continuation of the work begun by Sheldon and Ernst [14,15] in 
computer modelling of surface crack growth. The new version 
of the program attempts to overcome some of the problems that 
surfaced in the first program. There are also concepts 
integrated into the current program that were not in the 
original program. 



One new feature that has been introduced into the current 
version of the program is displacement control rather than 
load control. The current version of the program uses plastic 
displacement as the evolutionary variable. In the first 
program, the load was increased in small increments to 
simulate a mechanical test. One problem that arose from the 
use of this method was that virtually all the crack growth 
occurred in the last few load steps. Use of displacement 
control corrected this problem because with displacement 
control the load increases very quickly at first and very 
slowly around maximum load which is where the majority of 
crack growth occurs. Another feature of displacement control 
is that the simulation can continue past the maximum load; 
that is, the applied load can decrease. This feature would be 
useful in leak-before-break analyses. 

The program uses the key curve to calculate the load at 
each displacement step given the plastic displacement and the 
current crack geometry. The load and the crack geometry are 
then used to calculate the fracture variables along the crack 
front. 

The incrementing of plastic displacement is not intended 
to directly parallel an actual mechanical test. Plastic 
displacement was chosen because it was a nondecreasing 
variable in the key curve equation. The other possible 



evolutionary variables were load, the crack dimension a, and 
the crack dimension c. 

The program does not use any component of displacement in 
the calculation of the conditions along the front. Although, 
it is possible to calculate the elastic and plastic components 
of displacement in a real mechanical test, through unloading 
compliance measurements, in order to relate the simulation to 
an actual test. 

The original modelling technique used by Sheldon [15] 
divided the crack front into a collection of discrete points. 
Sheldon assigned values of a and c to each point along the 
front by fitting an ellipse to the point in question by using 
the position of the point, the position of the center of the 
ellipse, and the slope at that point as determined from the 
neighboring points. Sheldon then used these variables in a J 
calculation method as developed by McCabe, Ernst, and Newman 
[16] which was based on the Newman-Raju equations [7]. A 
value of J was calculated for each point along the front. A 
value of crack extension at each point was obtained from the 
material resistance curve, and the point was then advanced in 
the direction normal to the simulated crack front. 

The original program would divide the crack front into 
two halves when the program detected that the surface crack 
was mushrooming, that is, when the crack started growing in a 
non-self -similar way. The program considered the crack to be 



mushrooming when it detected a point with an approximately 
vertical slope. The section of the front from that point to 
the surface was modelled as a through crack with a crack 
length equal to the x coordinate of the point in question. 
The remainder of the crack front was modelled as a semi- 
elliptical crack as before. 

While this scheme showed promise, there were some 
complications that could be avoided by the use of a different 
model. The main problem was that calculating J using two 
different models led to problems in establishing continuity of 
J and dJ/da along the front. The program scaled the J values 
in order to maintain continuity of J, but no consideration was 
given to continuity of dJ/da. The discontinuity in dJ/da led 
to excessive crack growth at the point of the discontinuity. 

In part, the program developed in this chapter represents 
an evolutionary step forward. Modelling techniques that 
caused problems in the original program were modified. The 
use of two separate J estimation methods has been replaced by 
a single method that is used for all the points along the 
crack front. In addition to the evolutionary modifications, 
new features were added to the program to better model the 
fracture process. The primary addition to the program was the 
incorporation of constraint effects. The original version of 
the program as written by Sheldon [15] did not consider 



constraint effects. The techniques used to model crack growth 
will be described in the following sections. 


Computer Progr am Development 

This program models the evolution of the crack front in 
a displacement controlled test of a rectangular surface crack 
panel in tension. The program requires prior knowledge of the 
material deformation properties, the plane strain resistance 
curve, and the key curve for the specimen. The plane strain 
resistance curve is intended to represent the material 
fracture properties for the most constrained state of stress, 
and the key curve for the specimen will be used to relate the 
plastic displacement and current crack geometry to the load. 
The key curve was obtained using the same material and 
specimen geometry under blunt notch, or non-growing crack, 
conditions as explained in Chapter II. 

The program simulates a growing crack test by 
incrementing the plastic displacement. The user must enter 
the final plastic displacement along with the initial semi- 
elliptical crack geometry. The conditions along the crack 
front are recalculated and updated with each displacement 
step. The user is free to specify the number of displacement 
steps in any simulation. 


In this model, the crack front is represented by discrete 
points. The user is free to specify the number of points along 
the crack front. Specifying a large number of points allows 
the user to see how the front evolves in detail, but the 
program takes longer to run. The solutions calculated are 
insensitive to the number of points along the front; so, there 
is no computational advantage gained by specifying a very 
large number of points. 

It was sufficient to only consider one half of the front 
due to the symmetry of the problem. Boundary conditions on 
crack growth direction were applied at the surface and at the 
deepest point in order to enforce symmetry. Modelling only 
half the crack front halves the time it takes to run a 
simulation without affecting accuracy. 

This program was written in GW-^BASIC. A listing of the 
computer program is included in the Appendix. 

J Estimation along the Crack Front 

This program models the fracture process using a J 
resistance curve approach that incorporates constraint 
effects. This methodology was based on work done by McCabe, 
Ernst, and Newman [16]. In the first version of this program 
as developed by Sheldon and Ernst [14,15], a one-parameter 
approach to predicting fracture was used. In that program, J 
was calculated at each point along the crack front; then, the 



crack extension was calculated from the material R-curve. 
This approach did not consider constraint effects. In the 
current version of the program, a two-parameter approach 
incorporating constraint was used to calculate crack extension 
at each point along the crack front. The first step in 
modelling the crack extension across the crack front was to 
calculate the value of J at each point on the front. 

In Chapter I, it was shown that, given two cracked bodies 
subjected to the same loading history that differ in crack 
length by a small amount, the J integral can be defined as the 
difference in potential energy divided by the difference in 
cracked area. For a body where the crack only has a length 
dimension, the difference in cracked area, dA, is equal to the 
specimen thickness multiplied by the difference in crack 
lengths. It is not so simple for a geometry such as a surface 
crack that requires two length parameters. For example, a 
semi-elliptical surface crack requires the two length 
parameters a and c. 

For a semi-elliptical surface crack, Ernst [23] showed 
that the expression for the global value of J is linked to the 
way the virtual crack extension is taken. Ernst examined 
three different ways for taking the virtual crack extension: 
1) increasing c, keeping a constant, 2) increasing a, keeping 
c constant, and 3) increasing a and c, keeping a/c constant. 



The following equations for the plastic portion of J were 
developed by Ernst. The only difference between the three 
cases is the expression for the coefficient tj pl . 


J pi 




(3.1) 


Taking the key curve equation, which relates load and plastic 
displacement, and substituting it for the variable P yields 
the following expression. 


j pi 


= Pc 


IfiL 


( tc / 2 ) a c 


wcj layi layi ( p 

N+ 1 \ Cl \tj { WtpJ 


(3.2) 


The expressions for r^ pl turned out to be independent of the 
a/c and a/t ratios. These expressions are shown in Table 3.1. 


Table 3.1: Expressions for i) pl 


Variable 

i) plastic 

Held Constant 

Expression 

a 

-(p+m)/2 

c 

P/2 

a/c 

-in/ 2 




At this point it was necessary to determine which method 
should be used to calculate the global J and how that global 
j relates to the values of J along the front. The following 
normalization scheme was developed for this purpose. 

The total system energy in a cracked body can be 
expressed as follows. 

W = E + U (3.3) 

W = total energy supplied by external forces 
E = energy spent growing the crack 
U = strain energy 

The energy release rate for an incremental crack 
extension can be written as follows. 

EE = JL [ jds dn (3.4) 

bA bA J 

0 

The differential element ds represents an element along the 
crack front, and the differential element dn represents an 
element normal to the crack front. The above integral 
expression should be integrated along the entire crack front. 

In addition, taking the derivative of the energy balance 
with respect to the virtual increase in cracked area yields 
the following relationship. 



(3.5) 


6 W _ (6g + 5C7) 
6A " 6 A 


Replacing the incremental quantities in the above 
equation with experimentally measurable global quantities 
yields the following equations for constant displacement 
conditions . 


6A 



(3.6a) 


4 ^ = 4 - f Pdv (3.6b) 

6A dA J 
0 

The above result is entirely independent of the geometry of 
the surface crack. 

Now, for the case of growth of a semi-elliptical surface 
crack at a constant aspect ratio, Ernst shows that the 
following equations can be used to define ds,dn, and dA in 
terms of a, c, and elliptical angle, 0 [23]. 

ds - V(c cos<|>) 2 + (a sin4>) 2 dfy (3.7a) 


dn - 2 c JTc cos<|>) 2 + (a sin4>) 2 da 


(3.7b) 



dA = n c da 


(3.7c) 


Using these relationships allows the following definition 
of the energy release rate. 

£ 

|| - e * | / <#■ <3-8> 

This result shows that the global J obtained by taking 
the virtual crack extension with a/c constant is numerically 
equal to the linear average of J along the crack front. 

Since no assumptions were made regarding material 
behavior, the above equation applies for the elastic portion 
of J, the plastic portion of J, or the total J. 

This J estimation technique uses the observation that the 
distribution of the plastic portion of J along the crack front 
follows the distribution of the elastic portion of J for 
moderate amount on deformation [24,25]. This relationship can 
be expressed as follows. 

J pl (40 = k G(4>) (3.9) 


Therefore, 


J-(4>) = ( 1 + k) G(4>) 


(3.10) 



Following this line of reasoning the linear averages 


along the crack front should also be related as follows. 


k = 


J Jih±S. 




(3.11) 


Out of convenience, k was incorporated into a coefficient 
D that is independent of applied stress. 


D = 


k 

c (n ♦ 1) 


(3.12) 


By substitution, the value of D can be shown to be equal to 
the following expression. 


D = 


o . 224 n w 

* P 0 (n + 1) (G avs /o 2 ) \c) \ tj 


(3.13) 


The values of J along the crack front can be calculated using 
the following formula. 

J< 4») = G(4>) ( 1 + D a n - x ) (3.14) 

Now, the linear elastic portion of J as predicted by the 
Newman-Raju equations can be used with the above equation to 
perform an elastic-plastic J analysis along the crack front. 

Since the Newman-Raju equations are based on a semi- 
elliptical crack front, it was necessary to assign values of 



a/c and a/t to each point. A value of a/c was assigned to 
each point based on the a/c of an ellipse centered at the 
origin that passed through the deepest point and the position 
of the point in question. The value of a/t was taken to be 
the depth of the central point divided by the thickness. This 
method is illustrated in Figure 3.1. 

While this method of assigning values for a/c and a/t at 
each point is admittedly crude, the method was chosen because 
of its simplicity. The fact that this method is applicable at 
every point along the crack front prevents some of the 
problems that developed in the first attempt at simulating 
surface crack growth as developed by Sheldon [15]. 

Crack Growth Direction 

For the purposes of this computer program, crack growth 
was assumed to occur normal to the local crack front as 
defined by the points that represent the evolving crack front. 
The normal direction was determined by taking a numerical 
derivative at the point in question. 

As boundary conditions, the point at the front surface 
was constrained to growth along the surface, and the deepest 
point was constrained to growth through the depth. 



Ellipse fit to die three points 



Figure 3.1: Procedure for Fitting an Ellipse 


Crack Tip Constraint 

Today, it is accepted that constraint influences the J R- 
Curve. ASTM acknowledged this idea in their recommendation of 
bend specimens for determining R-curves [17]. This 
recommendation was based on the fact that bend specimens 
exhibit more in-plane constraint than tension specimens, and, 
thus, produce a more conservative design specification. It 
should be noted here that two specimens with identical out-of- 
plane constraint, e.g. plane strain or plane stress, can 
exhibit different levels of in-plane constraint. 



The typical differences in the R-Curves for bend and 
tension geometries illustrates how different levels of 
constraint affect the fracture properties of a material. The 
rate of change of the resistance curve slope is a function of 
constraint. Specimens under a high level of constraint have 
a slope that is lower than for low constraint configurations. 
Thus, the R-Curve for a high constraint situation would be 
below the R-Curve for a low constraint situation. 

It has been widely shown that a one parameter approach to 
modelling surface crack growth is inadequate [26,27,28]. The 
incorporation of a measure of constraint into this computer 
program is an attempt to address the problems inherent in a 
one-parameter approach. 

It has been shown by Ernst, Rush and McCabe [19] that 
there is a relationship between dJ/da and the level of 
constraint at the crack tip which they called (£ el ) -1 . The 
variable (£ e ^) 1 has been shown to characterize constraint. 
This variable was defined as follows. 



In calculating derivatives with respect to a virtual 
crack extension, a problem is again encountered in defining 
the manner that the virtual crack extension occurs. For the 
purpose of derivatives, the crack extension is assumed to 



occur normal to the equivalent ellipse that is fit at the 
point in question. This convention was adopted in the 
interest of consistency since the values for J are calculated 
from the equivalent ellipses. 

For the purposes of this computer simulation, dJ/da was 
chosen as a tool to incorporate constraint into the modelling 
of surface crack behavior. Briefly stated, the procedure 
calculates constraint effects by comparing dJ/da as calculated 
at a point on the crack front to dJ/da as calculated from the 
plane strain material resistance curve. 

Constraint effects were calculated in the following 
manner. The value of dJ/da was calculated for each point on 
the crack front. These values were obtained by calculating 
the value of J a small distance ahead of the point in 
question. Another equivalent ellipse was fit through the new 
incrementally advanced point. The a/c and a/t obtained from 
the three point fit were then used to calculate J. The value 
of dJ/da was found by dividing the difference in J values by 
the distance between the points. 

The methodology for incorporating constraint into the 
crack growth model began with the R-curve for (£ el ) _1 equal to 
zero. This R-curve represents the most constrained situation 
possible and served as the reference for the fracture 
parameters calculated in the program. The following equation 
was used to model the R-curve. 



A a r ./«. n « = «o cJ 3 ♦ J 2 + y 0 j 


(3.16) 


The R-curve for the side-grooved compact tension specimen 
was used in the computer program because the side-grooved 
compact tension specimen R-curve represented the most 
constrained geometry tested. The bend type loading and the 
side-grooves contributed to the high level of constraint in 
these specimens. 

The derivative of Equation 3.16 served a reference to be 
compared to the values of dJ/da calculated in the program. 


dJ\ m 1 

reference 3 tt c <7 2 + 2 0 o i7 + Y o 


(3.17) 


The factor F was defined as follows. 
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The prior assumption that constraint and dJ/da are 
related implies that constraint and the variable F are 
related. Substituting Equation 3.17 into Equation 3.18 yields 
the following. 



F 


(3.19) 


dJ 

da 3 <t 0 J 2 + 2 piJ + Y 0 


Substitution of the above equation into the integral 
definition for crack extension yields the following 
expression. 

Aa-/ SdJ (3.20) 

•o 


Integrating provides the formula for crack growth which 
incorporates constraint effects. 

Aa = A<3fa/afence (3.21) 

F 


Every point along the crack front had unique values for 
J and F. Together, the two parameters were used to calculate 
the crack extension at each point along the front. Once every 
point on the front had been advanced the appropriate amount, 
the plastic displacement was incremented, and the process 
began again at the new load. 



Comparison of Computer Output to Experimental Data 


It was initially hoped that this computer program could 
be used to model the 21-6-9 stainless steel surface crack 
tension specimens that were tested as part of this work, but, 
as shown in Chapter II, a J-based approach is invalid for this 
material. As a consequence, this computer program is 
inappropriate for prediction of crack growth in these 
stainless steel specimens. 

While this program is admittedly inappropriate for 
modeling the fracture behavior in these stainless steel 
surface crack specimens, it is convenient to use the material 
data developed in testing these specimens to test the general 
behavior of the program. The purpose for using the data from 
the stainless steel specimens is simply to check the program 
for instability, discontinuities, or any other obvious 
shortcomings . 

The results of three computer simulations are shown in 
Figure 3.2, 3.3, and 3.4. These figures compare the actual 
test results to the computer simulation for surface crack 
tension specimens SC16, SC11, and SC4, respectively. These 
surface crack specimens were chosen because they were all 
tested to a point past their maximum attainable load and they 
all exhibited significant crack growth. More specifically, 
these three specimens were chosen because their different 







Comparison of Simulation to Test Data 



Figure 3.2: Comparison of Simulation to Test Data for SCI 6 




Comparison of Simulation to Test Data 



Figure 3.3: Comparison of Simulation to Test Data for SC1 1 



Comparison of Simulation to Test Data 



Figure 3.4: Comparison of Simulation to Test Data for SC4 



aspect ratios. This allows an analysis of the performance of 
the computer program over a wide range of aspect ratios. 


Discussion and Conc lusions 

The results of the computer simulation using the data 
from the stainless steel surface crack tests shows that the 
program performs as expected. The new crack fronts calculated 
by the program do not exhibit any discontinuities or 
instability problems. In order to effectively evaluate the 
accuracy of the program, it would be necessary to compare the 
predicted results to the experimental results obtained from a 
set of surface crack specimens tested under J-dominant 
conditions. 

The experimental results presented in Chapter II show 
that the growth of these cracks is a deformation process. The 
photographs of these same specimens in Chapter II show the 
extensive plastic deformation around the crack mouth and at 
the back surface. The fact that the maximum loads for these 
specimens can be accurately predicted from deformation theory 
serves to reinforcing the assessment that a fracture mechanics 
based simulation for this material is inappropriate. 

It is not possible to judge the simulation results based 
on a comparison with the stainless steel surface crack 


specimens tested here. To make a definitive judgement on the 
accuracy of the simulation, it would be necessary to compare 
results for a material that was known to be under J-dominant 
testing conditions. 
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CHAPTER III: An Excerpt from WJ Curtin's Thesis 


CHAPTER HI 

FINITE ELEMENT INVESTIGATION OF 
CONSTRAINT PARAMETER 


Introduction 


Constraint is defined as the degree of crack tip stress triaxiality and is 
commonly quantified by the higher order crack tip stress and displacement 
terms (i.e. by T-stress or Q). Therefore, the elastic-plastic near tip fields are 
characterized by two parameters, J and a constraint parameter. The goal of the 
finite element analysis (FEA) investigation was to determine the coefficients of 
the higher order crack tip stress and displacement terms. As a subsequent step, 
the coefficients of the higher order terms were compared to the earlier proposed 
constraint parameter, (LJ 1 . The finite element investigation for this work is just 
the first step in a larger finite element investigation. That is, this work has 
started the development of an FEA data base containing the values of the higher 
order stress and displacement term coefficients. Future finite element 
investigation to determine these higher order coefficients shall include many 
other geometries. For example, three-dimensional finite element models of 
planar specimens will be developed and investigated. Upon completion of all 
FEA work the coefficients shall be used to compare various constraint 



parameters - including (L d ) ‘, q, Q, T-stress, and others — to the higher order 
terms. 

The coefficients of the higher order stress and displacement terms 
represent one part of a larger database which shall also include information on 
various constraint parameters and experimental results. The FEA results are an 
important part of the database; however, experimental results are critical to the 
evaluation of a constraint parameter. Experiments are very important because 
of the processes which lead to fracture. Ductile crack growth occurs due to 
micro void nucleation, growth, and coalescence. FEA cannot account for the 
processes which lead to ductile crack growth. Therefore, both experimental and 
FEA results are required in order to judge the ability of a second fracture 
parameter to characterize constraint. Each parameter shall be scrutinized upon 
completion of the entire database. 

Three types of FEA models were developed for this study -- a two- 
dimensional linear elastic center crack tension (CCT) model, a two dimensional 
elastic plastic CCT model and a three-dimensional linear elastic surface crack 
tension (SCT) model. Since this study has concentrated on elastic-plastic three 
dimensional cracks both a three dimensional model and an elastic plastic model 
were developed for the initial part of the finite element constraint parameter 
investigation. 



Finite Element Models 


Linear Elastic CCT Model 

All finite element work done for this study utilized ABAQUS software 
[28], The linear elastic CCT model consists of a two dimensional plane strain 
mesh using two planes of symmetry (Figure 3.1). That is, only one quarter of 


ThrMfh Cractd 



Figure 3.1: Center Crack Tension Model 


the specimen was modeled due to symmetry considerations. The mesh contains 
8-node biquadratic hybrid elements with reduced integration. The model 
consisted of 392 elements and was loaded in tension. A plot of the mesh and a 
sample of the ABAQUS input file are shown in appendix A* A total of five 
CCT meshes were developed for the analysis. The crack length over width ratio 
(a/w) varied from 0.25 to 0.75 for the five meshes (Figure 3.1). For consistency 
purposes, the same element arrangement near the crack tip was used for all five 
models. In order to correctly model the 1 / 4r stress singularity for the linear 


♦Appendix A not included 





elastic case, quarter point node elements [29] were used at the crack tip and all 
crack tip nodes were tied together. The element size (radial direction) at the 
crack tip was 0.0006e for the a/w=0.75 mesh where e=w-a. 

Linear Elastic SCT Model 

The linear elastic SCT model consisted of a three dimensional mesh 
using two planes of symmetry (Figure 3.2). The mesh contains 20-node 



Figure 3.2: Surface Crack Tension Model 

biquadratic displacement hybrid brick elements with reduced integration. The 
model consisted of 1764 elements and was loaded in tension. A plot of the 
mesh and a sample of the ABAQUS input file are shown in appendix B* In 
order to correctly model the 1 / 4r stress singularity for the linear elastic case, 
quarter point node elements were used at the crack tip and all crack tip nodes 
were tied together [29], The element size at the crack tip (radial direction) at 


Appendix B not included 



<J>=90° was 0.0024d where d=t-a. The modeled surface crack displays a shape of 
a/c=1.0 and a/t=0.5 (Figure 3.2). The modeled specimen thickness (t) was 1.0 
inch and the half width (w) was 2.0 inches. 

Elastic Plastic CCT Model 

The elastic plastic CCT mesh was identical to the linear elastic CCT 
mesh except for some adjustments. In order to model the r stress singularity for 
the elastic plastic analysis, quarter point node elements were not used at the 
crack tip and the crack tip nodes were not tied together. This allows for the 
crack tip to blunt for the elastic plastic case. The Ramberg-Osgood non-linear 
material deformation model was utilized for the elastic-plastic CCT model. In 
one dimension the Ramberg-Osgood model is, 

e a 
— = — + or 

*0 °"o 

where, 



and E is Young's modulus. The stress, strain, and yield strength are designated 
by ct, e, and ct 0 respectively. The material dependent constants a and n are the 
yield offset and hardening exponent. The material properties used in the model 
are displayed in table 3.1. These constants correspond to the properties of a 
typical low carbon steel. The element size (x direction) at the crack tip was 
O.OOle for the a/w=0.75 mesh where e=w-a (Figure 3.1). Similar to the linear 
elastic model, various a/w values were analyzed for the elastic plastic CCT 
case. Each model was loaded in three steps. The first load step went up to 70% 


' a ' 


(3.1) 



of limit load. The second step went up to 85% of limit load and the third step 
went up to 100% of limit load. A sample of an ABAQUS input file is shown in 
appendix C. 


Table 3 .1: Material Properties of a Typical Low Carbon Steel 


MATERIAL PROPERTY 

VALUE 

Young’s Modulus (E) 

30.0 E6 [psi] 

Yield Offset (a) 

0.5 

Strain Hardening Exponent (n) 

5 


40,000 [psi] 

Flow Stress (a f ) 

50,000 [psi] 


Theoretical Background 

In order to calculate the coefficients, the dependence on r must be 
determined for each of the higher order terms in the crack tip stress and 
displacement series (Figure 3.3). Williams [30] and Westergaard [2] have 
developed methodologies to determine the higher order terms dependence on r 
for linear elastic material properties. Sharma and Aravas [32] and Li and Wang 
[3 1] have formulated the higher order terms using an asymptotic analysis for a 
Ramberg-Osgood hardening material. 









Figure 3.3: The In-plane Crack Tip Stress and Displacement Locations 


Linear Elastic Background 

For the linear elastic case, the format of the higher order stress and 
displacement fields can be formulated using two early approaches developed 
by Williams [30] and Westergaard [2]. Williams derived expressions for the 
crack tip stress and displacement fields by the eigenfunction expansion method. 
Utilizing this method for a symmetric stress distribution (mode I) in an infinite 
body with traction free crack faces, the crack tip stresses at 0=0° and the crack 
tip displacements at 0=180° were derived. The results follow: 
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"'(r- 1800 ) = ^{[f^ + (r 4 ^} 3+ [f A+ (f- 4 ^} 5+ ''j (34) 

where T x , T y , b, and b 3 are constants that depend upon geometry and loading 
conditions. H is defined as v/(l+v). Poissons ratio is v and p is the shear 
modulus. Williams [30] showed that T„ and T y are zero for 0=0° provided that 
the crack faces are traction free. However, more recent investigators have 
determined that T x should not be zero [21]. Notice that the r w term in equations 
3.2 and 3.3 and the r M term in equation 3.4 both depend on b, and b y 

Westergaard [2] was the first to demonstrate that the crack tip stress 
fields can be derived for certain geometries by introducing an analytic function, 
Z(z), where z=x+iy and / = V-T . The normal stresses and displacements 
utilizing the Z function are. 


o a = ReZ-ylmZ' 
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where the bar over Z indicates integration with respect to z. Z is considered a 


function of the form. 
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where the c’s are complex numbers. 

The coefficients of the r w terms in equations 3.2 and 3.3 were defined as 


y K and y y , respectively. The coefficient of the r M term in equation 3.4 was 
defined as B. The relationship between the y's and B for a two dimensional 



3 / 


linear elastic stress state was determined using the Westergaard stress function. 
For mode I loading and y=0, the second terms on the right hand side of 
equations 3.5 and 3.6 drop out. Note that y x =y y when y=0. After expansion of 
the complex numbers, evaluation of the normal stress at 0=0°, and evaluation of 
the displacements at 0=180°, the relationship between y and B was determined 
for a two dimensional cracked body. This relationship was found to be. 


B=~A 


3 E' 


(3.9) 


where E-E for plane stress and E'=E/(l-v 3 ) for plane strain. For the linear 
elastic models a normalized B was defined as. 


= (3.10) 

and a normalized y as, 

r* = r^ (3.ii) 

where K is the stress intensity factor. By combining equations 3.9 to 3.11, the 
relationship between B* and y* was determined as. 


= (3.12) 

Rice [38] introduced the boundary layer approach to determine the crack 
tip stress and displacement fields. The boundary layer approach assumes that 
the boundary value stresses along the crack face are given by the extension of 
validity of the singular term in the elastic stress solution (Equation 1.1) to large 



values of r and small scale yielding. Larsson and Carlsson [21] introduced the 
modified boundary layer formulation for the crack tip stress fields. They 
included a constant T-stress term with the singular term in the boundary layer 
stress field formulation (T x in equation 3.2). For a two dimensional infinite 
body Larsson and Carlsson determined that only the o m T-stress term (TJ was 
non-zero. However, Nakamura and Parks [31] have determined that three T- 
stress terms are present in linear elastic three dimensional cracked bodies. They 
determined that the o n , o a , and x a crack tip stress fields required non-zero T- 
stress terms for regions behind the crack tip. For a two dimensional elastic- 
plastic cracked body, O'Dowd and Shih have demonstrated the existence of 
constant higher order terms in both the ct w and a n fields [22]. They introduced 
these higher order stress terms using the Q parameter. 

Elastic-Plastic Background 

For the elastic-plastic analysis, the format of the higher order stress and 
displacement fields have been formulated by Sharma and Aravas [32] and Li 
and Wang [33]. These researchers considered a two-dimensional crack problem 
in a homogenous Ramberg-Osgood elastic-plastic material. The first two terms 
of the crack tip stress series and the first three terms of the displacement series 
were formulated by assuming an asymptotic expansion for the stress terms in 
the form of, 



as r— >0 and s<t. Substituting equation 3.13 into the governing equations, the 
first two terms in the stress and the first three terms in the displacement fields 
were determined as. 
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where Q 0 is a dimensionless constant that controls the magnitude of the second 
terms, s = -l/(n+l), and e* 0) is an elastic strain resulting from the singular stress 
term. The governing equations include equilibrium, compatibility, and 
constitutive equations. The formulation of equations 3.14 and 3.15 neglect the 
effects of specimen geometry and far field loading and assume that the higher 
order terms are separable in r and 0. Using a Galerkin finite element technique 
t was determined to be 0.055 [32] for a two dimensional crack in plane strain 
with a material strain hardening exponent of 5. 


Procedure 


The field values of the normal crack tip stresses (c a and ct^) along r at 
0=0° and the crack mouth opening displacements (u y ) along r at 0=180° were 
calculated by each of the finite element models (Figure 3.3). The coefficients 



of each higher order term in the crack tip stress and displacement series were 
determined by curve fitting the FEA data. 

Linear Elastic Procedure 

The following formats which are similar to equations 3.2 through 3.4 
were utilized to curve fit the crack tip stress and displacement fields calculated 
by FEA for all linear elastic models: 

<r xx (r,0 0 )=D x r :i + T x + y z rl (3.16) 

o >y (r, 0° ) = D y r^ +T y + y y r* (3.17) 

1 3 

w y (r,180°) = Ar 1 + Br~ 2 (3.18) 
where r is the distance from the crack tip. The first terms (r l/J terms in 
equations 3.16 and 3.17 and r' n term in equation 3.18) were subtracted from the 
field values determined by FEA. The remaining higher order terms were curve 
fit using a least squares methodology. The curve fits yielded values for the 
coefficients T x , y x , T y , y y , and B. Note that the modified boundary layer 
approach was used to curve fit the higher order terms [21]. 

The range near the crack tip over which the stresses and displacements 
were curve fit needed to be determined. For the linear elastic models, the 
element stresses and nodal displacements (field values) very close to the crack 
tip which were affected by the quarter point node elements were omitted from 
the curve fitted data. For the linear elastic CCT model, the curve fitting range 



did not significantly affect the values of the coefficients. However, a criterion 
was used to determine the outer bound of the curve fitted data. The outer bound 
was determined using the T-stress. T-stress is the normal stress in the x- 
direction acting on the crack face at 0=180° [21]. The linear elastic CCT finite 
element analysis determined that o m was constant over almost the entire crack 
face. This T-stress value at 0=180° was then compared to the T x values from 
curve fits extending over various distances from the crack tip at 0=0°. The value 
of T x varied, although not significantly, depending upon the distance from the 
crack tip over which the higher order o a terms were curve fit. Therefore, the 
curve fit range which displayed the same T x value as the T-stress value along 
the crack face at 0=180° was used to determine the coefficients in equation 3.16 
for the linear elastic CCT model. The same curve fit range was then used for 
the higher order terms in equations 3.17 and 3.18. 

For the linear elastic SCT model, the in-plane T-stress at 0=180° was 
constant over much of the crack face. However, this T-stress did not coincide 
with the curve fitted T x values (0=0°) possibly because of the three dimensional 
stress state. That is, the out-of-plane stresses may have affected the T-stress 
values. Therefore, the outer bound was determined as the distance from the 
crack tip over which the curve fitted coefficients best represented the FEA field 
values close to the crack tip. The curve fitting outer bound was determined to 
be approximately 0.2a for all elliptical angles, where a is the surface crack 
depth. 

After the successful defense of this thesis, an error with the curve fitting 
format of equation 3.17 was discovered. The T y term should not be included in 



the curve fits since it would invalidate the traction boundary conditions along 
the crack face at 0 = 180°. The term was found to be zero for the linear elastic 
CCT model curve fits. However, T y was included in the curve fits for the linear 
elastic SCT model. The inclusion of T y is incorrect. See the Discussion section 
of Chapter III for further information pertaining to the incorrect format. 

Elastic-Plastic Procedure 

The following formats which are similar to the asymptotic solutions of 
equations 3.14 and 3.15 were utilized to curve fit the crack tip stress and 
displacement fields for the elastic-plastic models: 

cr a (r,0°) = D tpx r- y ^' ) +Q x r l (3.19) 

<jJr,0°) = D ipy r- ]l ^ + Q/ (3.20) 

u y (r, 1 80° )= + B, p r^' y "' + C„r' +1 (3.21) 

where r is the distance from the crack tip, s = -l/(n+l), and t = 0.055 for a two- 
dimensional crack in plane strain with n=5. Note that the strain hardening 
exponent does affect the exponent of r in all of the terms in equations 3.19 to 
3.21. The HRR terms (r 1 '"* 1 ) terms in equations 3.19 and 3.20 and r 1 '^ 1 ’ term in 
equation 3.21) were subtracted from the field values determined by FEA. The 
remaining higher order terms were curve fitted using a least squares 
methodology to determine the values of the coefficients Q x , Q y , B^, and C^. 



In developing a one or two parameter fracture methodology, one has to 
assume that the asymptotic solution on which the criterion is based provides and 
accurate description of the near tip stresses over distances that are sufficiently 
larger than the fracture process zone [32]. Hutchinson [34] suggests that one 
condition for a valid fracture methodology is, 

R>3S, 

where R is the radius of the zone of dominance of the methodology’s crack tip 
stress solution and is sufficiently larger than the fracture process zone. The 
crack tip opening is designated by 8 t . Shih [35] has shown that the crack tip 
opening can be determined by, 

S,=d,~ 

where d n is approximately equal to 0.5 for n=10 and 0.2 for n=3 [35]. By 
interpolation the approximate value of d n for n=5 is 0.3. By substitution of 8„ 
Hutchinson's condition becomes, 


R>3d,— 

Therefore, R extends a distance r/(J/o 0 ) » 1.0 from the crack tip when n=5. The 
radial distance from the crack tip is designated by r. The data inside R was 
omitted when curve fitting the elastic-plastic CCT stress and displacement data. 
The outer curve fit bound was r/(J/o 0 ) =15. The outer bound is sufficiently 
larger than R and is the approximate extent of the crack tip plastic zone for all 
models. 



Constraint Parameter Procedure 


In order to compare with the higher order term coefficients (LJ 1 was 
calculated at constant load by, 



AG 


G__Afl 


(G i+1 -G') 




(3.22) 


where G is the energy release rate and a is the crack length. Incremental steps 
are designated by i+1 and i. That is, a* 1 corresponds to an in-plane crack size 
slightly larger than a 1 . In order to attain values for G* 1 , all FEA models were 
modified so that the crack length was slightly larger than the original model. 
G, ve was defined as the average value of G* 1 and G or. 


G 


avt 


G" 1 + G' 
2 


(3.23) 


The values of G were calculated by ABAQUS [28] using the internal J-integral 
subroutine which calculates G using the domain integral method for all linear 
elastic models. For the elastic plastic CCT models G was calculated as J d using 
the EPRI [36] estimation scheme (see Equations 3.32-3.35). 




Results 


The coefficients were plotted as a function of a/w for the CCT models 
and elliptical angle (<J>) for the SCT model. In addition, the coefficients were 
compared with (LJ 1 . 

Linear Elastic CCT Model Coefficients 

G was calculated by ABAQUS using the domain integral method for 
each of the linear elastic finite element models. The ABAQUS value of G was 
within 1% of the handbook [24] value for all linear elastic CCT models. 

The crack tip displacement field values at 0=180° were determined using 
the linear elastic finite element models for a/w's ranging from 0.25 to 0.75. A 
sample of the crack tip displacements along with the singularity displacements 
(At 1 ' 2 term only) is shown in figure 3.4. After subtracting the singularity term in 
the displacement series, where for plane strain, 

A(9- 180°) = — J^(l- v) (3.24) 

H V 71 

the value of B was determined by curve fitting for each crack size, where jj. is 
the shear modulus. Note that B is the coefficient of the r w term in the 
displacement series defined by equation 3.18. B was then normalized as B* by 
use of equation 3.10. The results show an increasing trend with increasing a/w 
(Figure 3.5). The normalized coefficient values are shown in table 3.2. 
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Figure 3.4: Crack Tip Displacements for Linear Elastic CCT Model (a/w=0.25) 






Table 3.2: Higher Order Term Coefficients For Linear Elastic CCT Model 


a/w 

B* 

T x * 

y.* 

Yv* 


0.25 

-0.5275 

-0.6152 

1.396 

1.415 

2.328 

0.375 

-0.3482 

-0.3962 

0.939 

0.972 

1.894 

0.50 

-0.2615 

-0.2820 

0.671 

0.748 

1.832 

0.625 

-0.2157 

-0.2117 

0.462 

0.603 

1.951 

0.75 

-0.1737 

-0.1576 

0.079 

0.372 

2.520 


The crack tip normal stress values in both the x and y directions at 0=0° 
were also recorded from the linear elastic CCT finite element analysis. The 
finite element values along with the singularity values are shown in figure 3.6. 
The differences between the FEA and singularity values of the crack tip stresses 
are due to the low constraint CCT geometry. That is, the crack tip stresses 
deviate from the singular values for low constraint geometries. The singularity 
terms, D x r l/5 and D y r 1 ' 5 , where 

D x {6= 0°) = D y (0= 0°) = -JL (3.25) 

were subtracted from the FEA field values of and ct w , respectively. The 
form of equations 3.16 and 3.17 were used to curve fit the re mainin g higher 
order terms in order to determine T x , T y , y x , and y y values. y x and y y were 
normalized using equation 3.11. T x and T y were normalized using, 

l*=T ,^ f (3.26) 


1 





f 




where a w is the applied tensile stress. The curve fits determined that T y was 
zero for all crack lengths. This is consistent with the Larsson and Carlsson 
modified boundary layer approach for a two dimensional linear elastic stress 
state. That is, Larsson and Carlsson defined T y as zero in equation 2.19. The 
values of T x *, y*, and y y * are plotted versus crack length in figure 3.7 and their 
values are shown in table 3.2. 

Linear Elastic SCT Model Coefficients 

G, values for the SCT model were calculated along the crack front by 
ABAQUS. These values of G, compared well with the numerical equations for 
G, developed by Newman and Raju (Equation 1.6). Figure 3.8 compares these 
values of G,. 

The displacements and stresses calculated by FEA were transformed to 
the local in-plane coordinate systems at elliptical angles (<J>) ranging from 0° to 
90° (see Figure 1.4). The local in-plane coordinate system is shown in figure 
3.3 and is always perpendicular to the crack front. 

The crack tip displacement field values at 0=180° were determined using 
the three dimensional linear elastic finite element model for elliptical angles of 
0°, 22.5°, 45°, 67.5°, and 90°. An example of the FEA displacement field plotted 
versus distance from the crack tip is shown in figure 3.9. The first (singularity) 
term could not be subtracted from the FEA displacements because it is 
dependent on the unknown out-of-plane stress state (i.e. plane stress or plane 
strain). Along the front of a semi-elliptical surface crack the out-of-plane stress 
state varies from approximately a state of plane stress at the surface (<t>=0°) to an 
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Figure 3.7: Higher Order Stress Term Coefficients for the Linear Elastic CCT Model 
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Figure 3.8: Energy Release Rate for 3-D SCT Model 



mm 



approximate state of plane strain at <^90°. Therefore both A and B were 
determined at each elliptical angle by curve fitting where Note that A and B are 
displacement series coefficients defined by equation 3.18. The values of A 
determined by curve fitting were close to the plane stress value of A at <J>=0° and 
to the plane strain value of A at <f>=90°. The coefficient, B, was then normalized 
as B* by use of equation 3.10. The results show an increasing trend with 
increasing elliptical angle (Figure 3.10). The normalized coefficient values are 
shown in table 3.3. A possible source of error is introduced when determining 
the value of A by curve fitting. Error may be introduced because the A term is 
the dominant term in the two term series. Therefore, relatively small changes in 
A can introduce significant changes in B. 

The crack tip normal stress values in both the local x and y directions at 
0=0° were also recorded from the linear elastic SCT finite element analysis. An 
example of these stresses versus distance from the crack tip is shown in figure 
3. 1 1. The singularity terms, D x r w and D y r ,/I were subtracted from the FEA field 
values of a w and a n , respectively. Local values of the stress intensity factor, K„ 
were used to calculate the singularity terms of equation 3.25. The form of 
equations 3.16 and 3.18 were used to curve fit the remaining higher order terms 
in order to determine T„, T y , y x , and y y values. The coefficients, y x and y y , were 
normalized using equation 3.11 and T x and T y were normalized using equation 
3.26. The values of T x *, T y *, y*, and y y * are plotted versus elliptical angle in 
figure 3. 12 and are shown in table 3.3. Note that after the successful defense of 
this thesis the curve fitting format used to determine T y * and y* was found to 
be incorrect. See the Discussion at the end of Chapter ID for more information. 
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Figure 3.1 1 : Crack Tip Stresses for Linear Elastic SCT Model ( =45°) 
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Figure 3.12: Higher Order Stress Term Coefficients for the Linear 
Elastic SCT Model 




Table 3.3: Higher Order Term Coefficients for Linear Elastic SCT Model 


<L_ 

B* 

mm 

m 

Y.* 


mm 

0° 

-0.2672 

-1.1519 

-0.2457 

3.6601 

4.4192 

3.055 

22.5° 

-0.2175 

-0.3677 

0.1258 

.-0.8967 

1.0986 

2.621 

45° 

-0.1614 

-0.4080 

0.1700 

-0.5552 

1.3473 

2.275 

67.5° 

-0.1163 

-0.4295 

0.1849 

-0.7787 

1.1209 

2.014 

90° 

-0.0920 

-0.4360 

0.1640 

-1.0449 

1.2546 

1.910 


Accurate stress values at <J>=0° were not obtained from the SCT model because 
the finite element mesh was not fine enough in the z direction to account for 
surface effects. Therefore, the values of T* and y* were extrapolated between 
the elliptical angles of 0° and 22.5° in figure 3.12. In addition, a possible source 
of error was introduced during the evaluation of the singular stress terms. The 
local values of G, were calculated by ABAQUS. In order to determine K,, from 
G,, knowledge of the out-of-plane stress state is required, but is unknown. 
Therefore, the condition of plane stress was assumed during the evaluation of 
the singular stress terms for the SCT model. 

Elastic-Plastic CCT Model Coefficients 

All coefficients were determined for the elastic-plastic CCT models at 
70% of limit load (P 0 ). The total J was calculated by ABAQUS using the 
domain integral method for each model. The total J values were within 5% to 
10% of the EPRI estimation scheme [36] values of J for all models. The plastic 
























zone size determined by the equivalent stress being equal to ct 0 at 9=180° was 
observed to extend between r/(J /o 0 ) = 13 and r/(J/a 0 ) = 20 for all elastic-plastic 
models where r is the distance from the crack tip. 

The crack tip displacement field values at 9=180° were determined using 
the elastic-plastic finite element models for a/w's ranging from 0.25 to 0.75. An 
example of the crack tip displacements versus distance from the crack tip is 
shown in figure 3.13. After subtracting the HRR term (A^r 1 ""* 1 ' in Equation 3.21) 
from the FEA displacements, where 


d tp (0= 180°) = ae 0 


r J 

K as 0 o 0 I n 


n/(n+ 1) 

S,(l80°) 


(3.27) 


the value of and were determined by curve fitting. In equation 3.27, a 
and n are the yield offset and strain hardening exponent respectively. I n is an 
integration constant where I n = 5.0 for plane strain [37]. The material yield 

strength is designated by ct 0 and e 0 = aJE. J was calculated by ABAQUS and 
5,(180°) =2.3678 [37], and C^, were then normalized where, 
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(3.29) 


The coefficients, B^* and C^*, versus a/w are plotted in figure 3.14. The 
normalized values are shown in table 3.4. B„ and C were normalized by the 
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Figure 3.13: Crack Tip Displacements for Elastic-Plastic CCT Model (a/w=0.25) 







coefficients of the respective terms in equation 3.17 to obtain a dimensionless 
value. An additional l/(J/a 0 ) term was then included to obtain normalized 
units of [ 1/length]. The units of (LJ 1 are also [1/length]. 


Table 3.4 

Higher Ord 

er Term Coel 

ficients For 1 

Elastic-Plastic CCT Model 

a/w 

B„* 

C„* 

Q.* 

0* 

mm 

0.25 

176.62 

1206.2 

-1583.8 

-1391.9 

2.30 

0.375 

104.21 

1129.4 

-1284.7 

-1149.2 

1.90 

0.50 

61.83 

1265.5 

-1311.3 

-1213.5 

1.90 

0.625 

51.40 

1429.9 

-1447.5 

-1370.9 

2.13 

0.75 

46.81 

1793.5 

-1868.8 

-1809.7 

2.82 


The crack tip normal stress values in both the x and y directions at 0=0° 
were also recorded from the elastic-plastic CCT finite element analysis. An 
example of the normal stresses versus distance from the crack tip (r) is shown in 
figure 3.15. The HRR terms, D epjt r l,(inl) and D <py r l '°** l) , where 


A„(0=o°)= 


' j 



(3.30) 


were subtracted from the FEA field values of and c w , respectively. The 
values of the functions of theta were ct joc (o°) = 1.6836 and 5-^(0°) =2.2172 [37]. 


The form of equations 3. 19 and 3.20 were used to determine Q„ and Q y by curve 
fitting. All Q's were normalized in a similar fashion as using, 
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Figure 3.15: Crack Tip Stresses for Elastic-Plastic CCT Model (a/w=0.25) 


(3.31) 


Q* = 


Q, 



f-U 

W/cr 0 )) 


The values of Q x * and Q* are plotted versus crack length in figure 3.16 and are 
shown in table 3.4. Both Q x and Q y were found to be non-zero quantities and 
similar in magnitude. This is consistent with the findings of O'Dowd and Shih 
[22], These researchers introduced a Q parameter in order to extend J based 
fracture mechanics to include low constraint geometries. The Q parameter 
quantifies constraint by representing the second term in the near tip stress fields. 
Their investigation determined that the higher order stress terms (Q x r and Q y r 
in equations 3.21 and 3.22) were approximately constant but did show a slight 
dependence on r for stresses in the forward sector of the plastic zone (-90° < 0 < 
90°). The present findings also demonstrated that the higher order stress terms 
showed only a slight dependence on r. 


Constraint Parameter Comparison 

The coefficients of the higher order crack tip stress and displacement 
terms were plotted with (LJ- 1 for each of the finite element models. These 
figures were developed as one part of a growing database of information on 
crack tip constraint. The database also includes information on additional 
constraint parameters and experimental results of various fracture geometries. 
The ability of these parameters to characterize constraint shall be scrutinized 
upon the completion of the entire database. 

(LJ 1 was calculated for each of the linear elastic CCT models. B* was 
then compared to (LJ 1 . Figure 3.17 and 3.18 show that -B* and -T* display 
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Figure 3.16: Higher Order Stress Term Coefficients for Elastic-Plastic CCT Model 
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Figure 3.17: (Lel)-I and -B* for Linear Elastic CCT Model 
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Figure 3.18: (Lel)-I and -Tx* for Linear Elastic CCT Model 



similar trends when plotted with (LJ 1 from a/w=0.25 to a/w=0.5, but at higher 
values of a/w the trends diverge. Figure 3.19 displays y* and y y * with (LJ 1 . 

For the linear elastic SCT model, local values of the proposed constraint 
parameter were calculated by incrementing the size of the crack and applying 
equation 3.22. Local values of virtual crack extension perpendicular to the 
original crack front were used to calculate (L d )'. Figure 3.20 displays similar 
trends for -B* and (LJ '. Accurate stress values at <|>=0 # were not obtained from 
the SCT model because the finite element mesh was not fine enough in the z 
direction to account for surface effects. Therefore, the values of T* and y* 
were extrapolated between the elliptical angles of 0° and 22.5° in figures 3.21 
and 3.22. T x * and -T y * show similar trends with the proposed constraint 
parameter when plotted as a function of elliptical angle (Figure 3.21). Figure 
3.22 compares y* and y* with (LJ ', but a similar trend was not found. Note 
that an error was found with the higher order stress term coefficients for the 
SCT model which is explained in the Discussion of Chapter III. 

For the elastic-plastic models, the linear elastic portion of J (G or JJ was 
calculated using the EPRI [36] estimation scheme, where 


J,= 


(1 -S)K’(aJ 


(3.32) 


a 'g =a + 


and, 

1 1 fn-lY*(a)Y 


1 + {P!PJ P* 
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(3.33) 
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Figure 3.19: (Lel)-I and Gamma* for Linear Elastic CCT Model 
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Figure 3.20: (Lel)-I and B* for Linear Elastic SCT Model 
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Figure 3.22: (Lel)-I and Gamma* for Linear Elastic SCT Model 



(3.35) 


where |3=6 for plane strain, P is the total applied load, and n is the strain 
hardening exponent. The specimen thickness, length of the remaining ligament, 
ultimate tensile strength, and yield strength are designated by t, b, a u , and cr 0 
respectively. Loaded to 70% of limit load (P # ), the ABAQUS calculated total J 
values deviated from the EPRI J d values by 15% to 35%. Therefore, 15% to 
35% of the total J was due to plasticity (J^. For the elastic-plastic CCT model, 
B^* and C^* were then plotted with (LJ 1 in figure 3.23 and 3.24, respectively. 
Figure 3.25 displays Q x * and Q y * with the proposed constraint parameter. 


Discussion 


The coefficients of the higher order terms developed from the linear 
elastic CCT model were compared with higher order term relationships derived 
using the eigenfunction expansion method [3,30] and the Westergaard [3] Z 
function (Equation 3.12). As predicted by the eigenfunction expansion and the 
Z function methods, y * is virtually equal to y y * in our numerical results. 
However, the values begin to deviate at higher a/w’s; additional higher order 
stress terms are possibly required with larger cracks since the finite body effects 
(at the surface) become more significant. The relationship determined by the 
Westergaard approach (Equation 3.12) between y* and B* was: 
B*=(-l/3)y x *=(-l/3)Y y , '‘. However, the FEA results (with a/w<0.5) indicate a 
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relationship of the form: B* * -0.37 y z * * -0.37 y y *. The difference in the two 
results could be attributed to finite body effects. 

After the successful defense of this thesis, the curve fitting format of 
equation 3.17 was found to be incorrect. The T y term was included when curve 
fitting the higher order terms for the linear elastic SCT model. However, the 
T y term should not have been included in the curve fits since it would invalidate 
the traction boundary conditions along the crack face. The term was not found 
to be zero, as it was in the linear elastic CCT model, because the higher order 
Gyy terms converged to a finite stress value at the crack tip as r— >0. The most 
likely explanation for the finite value convergence is the inaccuracy of the 
singular term calculation. An assumption on the out-of-plane stress state 
was made during the calculation of the singular term in order to convert the 
energy release rate, G„ calculated by ABAQUS [28] to the stress intensity 
factor, K,; after converting, K, was then used to calculate the singular stress 
terms. Perhaps a better procedure to determine the singular stress term includes 
iteratively offsetting the higher order terms by T y . That is, determine the 
higher order term coefficients using the out-of-plane stress state assumption. 
Then, add T y to the singular stress term to determine the value of an iterative 
singular stress term. At this point, the iterative singular stress term is subtracted 
from the FEA stresses to obtain iterative values of the higher order o n stress 
terms. The iterative higher order terms can be curve fit in order to obtain 
iterative values for T y and y y . Repeat the previous steps until the iterative 
singular a n term converges and T y =0. After convergence, a more accurate value 
for K, can be calculated using the singular stress term. This value of K, should 



then be used to calculate the singular term and the first term in die 
displacement series, so that the higher order and u y term coefficients can be 
determined. In addition, perhaps a similar iterative procedure could be utilized 
to bring the higher order term coefficients of the linear elastic CCT model 
(a/w=0.625,0.75) into agreement with equation 3.12. 

The Q values for the elastic-plastic CCT model were determined by 
curve fitting the two-term stress series expansion (Equations 3.19 and 3.20). 
The two terms yielded good representations of the crack tip stresses over the 
curve fitted range 1 < r/(J/a 0 ) < 15. Sharma and Aravas [32] used this two term 
series in their asymptotic analysis of crack tip fields. However, they also stated 
that at the distances r/(J/cr 0 )=2 and r/(J/crJ=5 from the crack tip additional 
terms in the asymptotic stress expansion (Equation 3.13) may be needed for an 
accurate representation of the stress field in front of the crack tip. 

The second and third order terms were required for curve fitting the 
elastic-plastic displacements. The Sharma and Aravas and Li and Wang 
analyses both state that the second order term (B^r*"- 1 ** 1 ) does not account for 
the effects of elasticity. The third order displacement term is a function of the 
elastic strain resulting from the singular (HRR) stress term. Therefore, with the 
second and third order terms included, excellent agreement between the FEA 
field values and the curve fitting constants were obtained for the displacement 
fields. 

The similar trends observed between C^* and (LJ-' in figure 3.24 and Q* 
and (LJ 1 in figure 3.25 are significant developments. The similar trend 
suggests that C^*, Q*, and (LJ 1 are related. Of course, an implicit relationship 



exists between the terms in the crack tip stress, strain and displacement series. 
Therefore, if a relationship does exist between (LJ ' and the higher order terms, 
(L el ) 1 can be considered a second fracture parameter (i.e. constraint parameter). 
However, an extensive analytical study of the higher order term coefficients is 
required in order to develop a relationship with (LJ 

(LJ 1 is a parameter associated with the higher order strain/displacement 
terms for a linear elastic stress state. C * is a function of the elastic 

«P 

strain/displacement resulting from the HRR term in the stress series. Therefore, 
the similar trend between C^* and (LJ 1 suggests that constraint parameters may 
be separable in elastic and plastic components by relating these separate terms 
to different higher order terms in the stress, strain, or displacement series. 
However, finite geometry effects have not been accounted for in the comparison 
between and (LJ 1 . In addition, experimental results must be implemented 
in the investigation before concluding that (LJ 1 is a valid constraint parameter 
for general use. Nevertheless, the possible relationship with the higher order 
stress and displacement terms suggests a bright future for (LJ 1 as a constraint 
parameter. 
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EFFORTS IQ CHARACTERIZE THREE-DIMENSIONAL EFFECTS 
OBSERVABLE HJ TS2Q DIMENSIONAL FRACTURE SPECIMENS. 

by D. Lambert and H. Ernst 
4.1. Introduction 

The underlying purpose of this research is to develop a 
methodology which would allow the characterization of three- 
dimensional (3D) effects in fracture. This characterization 
should include: 

(1) geometric effects arising from crack front 
curvatures (curvatures are present, for example, 
with surface cracks) , 

(2) geometric effects related to thickness and ligament 
length (the gross sizing details that affect the 
three-dimensionality of the state of stress at the 
crack tip) , and 

(3) loading geometry effects (including three- 
dimensionality of the far field stress arising from 
the character of applied loads, and, especially the 
gradient of the far-field stress arising from 
differing ratios of bending-to-tension) . 

Two-dimensional (2D) or planar specimens have been observed 
to generate different fracture resistance curves when different 


4.2 


thicknesses are tested. Specifically discussed here are J M R 
curves that use the J-modified parameter as developed by Ernst 
[ 4 . 1 , 4 . 2 ] ^ . Different configurations have been shown to support 
a differing degree of triaxiality of the stress field in the 
vicinity of the crack front, where the fracture process is 
occurring. The degree of stress field triaxiality that is 
exhibited is referred to as the constraint. The J^R curves are 
a result of the different, averaged constraint in each specimen. 
Even though the configurations are considered to be planar, 
curvatures can develop in crack fronts that result from fracture 
in the presence of a gradient of the constraint within the 
specimen. Thus, the complexities that occur in the most general 
cases of fracture appear in the simplest cases of planar 
specimens. Ultimately, to evaluate fracture resistance retaining 
a planar analogy requires that the crack front fall within 
specific limits of straightness. 

Since, 3D stress fields are present in planar 
configurations, an effort to map the crack face separation 
profiles of a variety of geometries as a function of the position 
within the cross-section has been proposed. The crack tip 
opening displacement (CTOD) is a linear function of the J- 
integral [4.3]. Profiling represents an extension of that 
functional relationship. 

One goal of the overall research program is to test a wide 
variety of planar specimens, varying the thickness and length of 
the initial remaining ligament, as well as the bending-to-tension 
•Numbers in brackets refer to the references at the end of the paper 
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ratio due to the nature of the applied load to produce 
significant changes in fracture behavior. The results will be 
compared using the ligament dimensions as variables, and finite- 
element analysis will be used to evaluate the triaxiality of the 
stress field. Although other parameters might be used, one 
parameter being considered to quantify the stress field 
triaxiality is h: 

fa _ Oman ( 4 . 1 ) 

OvM 


Here, a mean is the mean stress and a vM is von Mises equivalent 
stress. Ultimately, this approach is expected to produce 
parameters and fracture behavior that can be generalized to the 
3D cases that are of the most interest. 


4.2 introduction to Profi ling 

The displacement of the faces of a crack are a function of 
the loading and of the position. Using a two dimensional 
analogy, the displacement of a point along a crack face within an 
elastic body was given by Tada, et al [4.4] : 


v 


el 


E'yfn 


Kyfr 


(4.2) 
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Here, v e i is the "y" -directed elastic displacement at position r, 
measured from the crack tip to the point in question, E' is the 
equivalent modulus (E r =E for plane stress, E =E/ (1-v^) for plane 
strain, v is Poisson's ratio) . The loading is specified in the 
presence of a flaw by K, the stress intensity parameter. This 
displacement relationship has a square-root of r form. 

Hutchinson [4.5] and Rice and Rosengren [4.6] developed a 
similar form for plastic response that follows Ramberg-Osgood 
deformation characteristics, i.e.: 

— = «(—)" (4.3) 


In this equation, E and c are the equivalent strain and stress 
and E 0 , o 0 an< 3 n are mat erial constants. The form of the 
displacement is as follows: 


v pl =k-J%"r } 


(4.4) 


This equation is written for a non-growing crack, and v p i is the 
the displacement of the body, assuming Ramberg-Osgood type 
deformation, and k includes the functionality with regards to the 
constraint, i.e., plane stress or plane strain condition. 

Looking at the equations (4.2) and (4.4), the constraint appears 
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in the coefficient E' for the linear elastic case and in the 
coefficient k for the plastic case. Thus, the separation at 
various points through the ligament thickness could be expected 
to reflect that difference in constraint that arises with the 
depth into the thickness. It may also provide a measure of that 
constraint . 

Since the development above is for a non-growing crack 
situation, differences that occur between the theoretical elastic 
plus plastic displacements and the displacement profile of an 
actual growing crack near the crack tip might provide a fracture 
criterion on that local level. 

4.3 Profiling Matrix and Petflils 

One objective in the research was to characterize the 
separation between the surfaces of cracks. This separation 
profile is a function of the level of J and of the position 
within the ligament. In this case the position would include the 
distance from the load-line, x, in the direction of crack growth, 
and the depth beneath the surface in the thickness direction, z. 

The primary effort in the past six months has been to 
characterize the crack face separation of selected specimens. 

The specimen identities and the corresponding configurations 
appear in Table 4.1, below. Data was generated for a total of 
six compact tension (CT) specimens and three center-crack tension 
(CCT) specimens. The CCT configuration produces two crack 
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Table 4.1: Matrix of Specimens Profiled 


(W, B, and b are in inches, a/W is nondimens ional ) 


Spec 

Conf ig 

w 

B 

b 

a/W 

Remarks 

No of 
Profs 

E2 

1T-CT 

2 

1/2 

1 

. 5 

Baseline Specimen 

8 

81 

1T-CT 

2 

1/2 

.5 

.75 

Larger Init. a/W 

8 

51 

1T-CT 

2 

1/4 

1 

.5 

Thinner Section 

5 

C9 

1T-CT 

2 

.85 

1 

. 5 

Thicker Section 

13 

82 

1T-CT 

2 

1/2 

1 

.5 

Multi -Specimen 

8 

84 

1T-CT 

2 

1/2 

1 

. 5 

Multi -Specimen 

8 

08 

1T-CT 

20%SG 

2 

1/2 

.8 

.6 

Side-Grooved 

8 

D6 

CCT 

1 

1/2 

.5 

. 5 

Tension 

2x8 

55 

CCT 

1 

1/8 

.5 

. 5 

Tension, Thin 

2x3 

B9 

CCT 

1 

.85 

.5 

.5 

Tension, Thick 

2x9 


profiles per specimen, and thus the total number of crack fronts 
observed is thirteen. Profiles were made of the AL6061-T651 in 
every case: the IN718-STA1 material has proven too hard to 

polish in the same fashion as the aluminum. Epoxy infused into 
the gap of the crack was effective in producing a well-defined 
crack profile for the aluminum, but the profiles of the nickel 
were rounded and poorly-defined. Until the techniques have been 
modified to overcome the rounding, the profiling was suspended 
for the nickel. 

After mounting the specimens in epoxy, the exposed surface 
was polished to provide a surface profile. After recording the 
profile, 0.025- to 0.035-inches was removed by grinding and 
polishing to produce the next profile to be recorded. This was 
continued for each specimen into the center of the cross-section 
























































































14440 REM 

15000 REM dJ/da CALCULATION 

15001 REM 

15100 REM TRANSFER ARRAYS SO EXISTING SUBROUTINES CAN BE USED 
15105 REM 

15110 FOR DUM-0 TO NUM: JT(DUM) -J(DUM) : XT(DUM) -X(DUM) : YT (DUM) =Y (DUM) : NXT (D' 
DUM) :NYT(DUM)-NY(DUM) :NEXT DUM 

15120 REM 

15150 REM CALCULATE NORMAL TO EQUIVALENT ELLIPSE 
15155 FOR DUM=1 TO NUM-1 

: 15160 YD (DUM) — X (DUM) *A* ( 1- (X (DUM) /C(DUM) ) * 2 ) “- . 5/C ( DUM) *2 
15162 PX— -YD (DUM) :PY-1 
15164 HYP— ( 1+PX * 2 ) * . 5 
15166 NX (DUM) -PX/HYP 
15168 NY (DUM) -PY/ HYP 
15195 NEXT DUM 

15200 REM CALCULATE COORDINATES OF POINTS AHEAD OF FRONT 

15210 FOR DUM-0 TO NUM 

15220 X (DUM) -XT (DUM) +NX(DUM) * (T/2500) 

15230 Y ( DUM ) - YT ( DUM ) +NY (DUM) * (T/2500) 

15240 NEXT DUM 

15241 A-Y (NUM) 

15242 REM 

15245 REM CALCULATION OF EFFECTIVE C'S SUCH THAT NORMALS ARE EQUAL 
15250 FOR DUM-1 TO NUM-1 

15255 C (DUM) -(X (DUM) *2/ ( 1- (Y (DUM) /A) *2) )'.5 
15260 NEXT DUM 

15265 C (NUM) =C (NUM-1) :REM BECAUSE C AT THE TOP IS UNDEFINED 
15270 C ( 0) -X ( 0) : REM BY DEFINITION 

15301 REM 

15302 REM CALCULATE ELLIPTICAL ANGLES 

15304 GOSUB 12000 

15310 REM 

15320 REM CALCULATE G'S 

15330 GOSUB 8000 

15340 REM 

15350 REM CALCULATE J'S - 

15355 GOSUB 8500 

15360 REM 

15363 REM CALCULATE dJ/da FROM INCREMENTAL GROWTH OF EQUIV. FRONTS 

15364 REM 

15365 FOR DUM-0 TO NUM: JDA (DUM) = (JT ( DUM) -J (DUM) )/ (T/2500) : NEXT DUM 

15390 REM 

15400 REM CALCULATE dJ/da FROM MASTER R-CURVE 

15410 FOR DUM-0 TO NUM 

15420 JDAM(DUM) — 3*ALPA*J (DUM) *2+2*BETA*J (DUM) +GMMA 
15425 JDAM (DUM) — 1/ JDAM (DUM) 

15430 NEXT DUM 

15450 FOR DUM-0 TO NUM 

15455 F( DUM) -JDA (DUM) /JDAM (DUM) 

15456 IF F (DUM) <1 THEN F(DUM)«1:REM BY DEFINITION F CANNOT BE LESS THAN 1 
15460 NEXT DUM 

15500 REM TRANSFER BACK ARRAYS 

15510 FOR DUM-0 TO NUM : J (DUM) -JT (DUM) : Y (DUM) -YT (DUM) : X (DUM) -XT ( DUM) : NX (DUM 
DUM) : NY (DUM) -NYT (DUM) : NEXT DUM 
15600 RETURN 



10020 

10040 

10050 

10200 

10210 

10220 

10900 

10998 

10999 

11000 
11001 
11002 
11100 
11150 
11200 
11250 
11300 
11495 

11499 

11500 
11550 
11600 
11650 
11700 
11900 

11997 

11998 

11999 

12000 
12001 
12100 
12110 
12115 
12117 
12120 
12150 
12400 
12450 

12998 

12999 

13000 

13001 
13100 
13110 
13120 
13130 
13140 
13150 
13160 
13400 
13450 

14000 

14001 
14100 
14105 
14110 
14120 
14200 
14420 
14430 


C (DUM) * ( (X(DUM) "2*A"2) / (A"2-Y(DUM) *2) ) ".5 
NEXT DUM 
C (NUM) *C (NUM-1) 

FOR DUM=0 TO NUM 

IF C (DUM) >W/ 2 THEN C (DUM) *W/ 2 : PRINT” C > W/2 
NEXT DUM 

R ™Ti*tt**t*t***tt*************************************A****** 


REM 


RE M CALCULATION OF Fs AT ONE POINT 


REM 

REM FOR A/Ol 

Ml* (C (DUM) / A) A .5*(1+.04*C(DUM) / A) 

M2*. 2* (C(DUM) /A) *4 ^ 

LG~1+ ( ^1+^35* (C (DUM) /A)* (A/T) " 2 ) * ( 1-SIN (TH (DUM) ) > * * 

FT=( (C(DUM) /A) * 2* (SIN (TH (DUM) } ) * 2+ (COS (TH (DUM) ) ) 2) .25 

RETURN 

REM 

Ml=1.13“9. 000001E-02* (C ( DUM) /A) 

M2=-.54+.89/(.2+A/C(DUM) ) 

M3= . 5-1 / ( .65+A/C(DUM) ) +14 * ( 1-A/C (DUM) ) "24 
LG=1+ ( . 1+. 35* (A/T) "2 ) * (1-SIN (TH(DUM) ) ) "2 
FT=( ( A/ C ( DUM) ) "2*COS (TH (DUM) ) "2+SIN (TH (DUM) ) "2) 


,25 


"S.V.nuntiiuHiiiiiiiintttinntimnHuntuHnnnn 


REM 

REM 


CALCULATION ELLIPTICAL ANGLE ALONG THE FRONT 
REM ~ 

IF A/C ( DUM) <1 OR A/C (DUM) =1 THEN OPP=Y (DUM) :ADJ* (A~2-OPP 2) .5 

IF A/C (DUM) > I THEN ADJ=X (DUM) : OPP= (C (DUM) * 2-ADJ'2 ) * . 5 

IF A/C ( DUM) <0 OR X (DUM) >C(DUM) THEN PRINT"A/C OUT OF RANGE :GOTO 500 
TH ( DUM ) =ATN ( OPP / ADJ ) 

NEXT DUM 


REM 

REM 


REM CALCULATION OF Dstar 

REM 

FOR DUM=0 TO NUM 
XD1=.224*N*W 

XD2*PI*BN* (N+l) * (GAVG/SIG*2) 

XD3= ( A/C ( DUM) ) * (1-PP*N) 

XD4* (A/T) * (-2-M*N) 

XDS ( DUM) =XD1*XD3*XD4 /XD2 
NEXT DUM 
RETURN 
REM 
REM 

REM INCREMENT PLASTIC DISPLACEMENT 

VPL-VPL+IVPL 

IF VPL>-VPLM THEN 500 

P-W*T*BN* (Y (NUM) /X(0) ) *PP* (Y (NUM) /T) ~M* (VPL/T) 
SIG-P/ (W*T) 

RETURN 

JDAM ( DUM ) = 1 / JDAM ( DUM ) 

NEXT DUM 


UHtUt*************************************************** 


(1/N) 



8498 REM 

8499 REM 

8500 REM POINT J CALCULATION 

8501 REM 

8502 REM 

8580 REM FIND GAVG 

8590 GTOT=0 

8600 FOR DUM=0 TO NUM 

8610 GTOT=GTOT+G (DUM) 

8620 NEXT DUM 

8630 GAVG=GTOT/ (NUM+1) 

8640 REM 

8700 GOSUB 13000 
8705 REM 

8710 FOR DUM=0 TO NUM 

8720 J (DUM) *G (DUM) * ( 1+XDS (DUM) * (SIG/BN) * (N-l) ) 

8725 NEXT DUM 

8995 REM 

8996 RETURN 


8997 REM ******** ft 

8998 REM 

8999 REM 

9000 REM DELTA-a CALCULATION 

9001 REM 

9002 REM 

9010 FOR DUM = 0 TO NUM 
9015 START ( DUM ) * ADEL ( DUM ) 

9020 ADEL (DUM) « (ALPA* J (DUM) A 3+BETA* J (DUM) A 2+GMMA*J (DUM) ) /F(DUM) 

9025 IF ADEL ( DUM )< START (DUM) THEN ADEL ( DUM )*= START ( DUM) 

9026 INCA ( DUM ) =ADEL ( DUM )- START (DUM) 

9027 REM 
9090 NEXT DUM 

9100 REM MOVING AVERAGE TO SMOOTH OUT INCREMENTAL CRACK EXTENSION 
9110 REM 

9170 FOR DUM=1 TO NUM-1 

9180 INCA (DUM) = ( INCA ( DUM- 1) +INCA(DUM) +INCA(DUM+1) ) /3 
9190 NEXT DUM 

92 0 0 INCA ( 0 ) = ( 2 * INCA ( 0 ) + INCA ( 1 ) )/3 

9210 INCA (NUM) = ( 2 * INCA (NUM) +INCA (NUM-1) ) / 3 

9496 RETURN 

9497 REM *********************************** ************************ 

9498 REM 

9499 REM 

9500 REM ADVANCE THE CRACK FRONT 

9501 REM 

9502 REM 

9510 FOR DUM=0 TO NUM 

9520 X (DUM) =XO (DUM) +NX (DUM) *INCA(DUM) : XO (DUM) =X ( DUM) 

9530 Y (DUM) =YO (DUM) +NY ( DUM) *INCA(DUM) : YO (DUM) *Y (DUM) 

9535 IF X (DUM) >W/ 2 OR Y(DUM)>T THEN PRINT: PRINT 11 EXCESSIVE CRACK GROWTH" 
9540 NEXT DUM 
9590 REM 

9997 RETURN 

9998 REM 

9999 REM 

10000 REM CALCULATION OF EFFECTIVE C ALONG CRACK FRONT 

10001 REM 

10002 REM 
10004 A=Y (NUM) 

10010 FOR DUM= 0 TO NUM-1 



4295 REM Httttttt**tttt**tt*tt*tt***t*********t******************** 

5998 REM 

5999 REM 

6000 REM PLOT THE CRACK FRONT ON SCREEN 

6001 REM 

6002 REM 

6100 LOCATE 1 , 1 : PRINT : PRINT "LOAD - ";P;" lbf Vpl - ";VPL;" in 

tl 


6110 

6125 

6130 

6140 

6150 

6996 

6997 

6998 

6999 

7000 

7001 

7002 

7003 

7004 

7005 

7006 

7007 

7008 
7010 
7020 
7025 

7029 

7030 

7035 

7036 
7040 
7045 
7050 
7060 
7065 
7070 

7996 

7997 

7998 

7999 

8000 
8001 
8002 
8010 
8020 
8025 
8330 
8335 
8400 
8405 
8410 
8430 
8440 
8450 
8460 

8496 

8497 


REM 

LINE (2.6*X(0) /T,Y(0) /T) - ( 2 . 6*X ( 1 ) /T, Y(l) /T) 

FOR DUM = 1 TO NUM- 1 

LINE ( 2 . 6*X (DUM) /T , Y (DUM) /T) - (2 . 6*X(DUM+1) /T , Y (DUM+1) /T) 

NEXT DUM 
RETURN 

REM *********************************************************** 


REM 

REM 

REM 

REM 

REM 


NY ( 0 ) =0 : REM 
NX ( 0) =1 : REM 
REM 

NX (NUM) =0: REM 
NY (NUM) =1: REM 
REM 


DIRECTION OF CRACK EXTENSION 


CONSTRAINS EDGE TO HORIZONTAL FREEDOM 
CONSTRAINS EDGE TO HORIZONTAL FREEDOM 

CONSTRAINS CENTER TO VERTICAL FREEDOM 
CONSTRAINS CENTER TO VERTICAL FREEDOM 


FOR DUM=1 TO NUM-1 
IY«Y(DUM+1) -Y (DUM-1) : REM 
IX=X (DUM+1 ) -X (DUM-1) : REM 
REM 

PX=IY :REM 
PY«- IX:REM 


DIFFERENCES IN POSITIONS FOR NEIGHBORING 
POINTS ON THE CRACK FRONT 

CONVERT TO MEASURE OF SLOPE OF 
THE PERPENDICULAR 


REM 

HYP=(PX*2+PY*2) * . 5 
REM 

NX (DUM) =PX/HYP 
NY (DUM) =PY/HYP 
NEXT DUM 


REM 

RETURN 

REM *********************************************************** 


REM 

REM 

REM 

REM 

REM 

FOR DUM=0 TO NUM 

PHI 2® ( 1+1 . 464* (C ( DUM) /A) "1.65) : REM 
REM 

FW= ( COS ( (PI*C (DUM) /W) * (A/T) ".5) ) "-.5: REM 
REM 

IF A/C (DUM) >1 THEN GOSUB 11000 ELSE GOSUB 11500:REM Ml , M2 , M3 , LG , FT 
REM 

FS=(M1+M2MA/T) *2+M3* (A/T) *4)*LG*FT*FW 
REM 

K«SIG*FS*(PI*A/PHI2) *.5 
G (DUM) =K* 2 / E 
NEXT DUM 
RETURN 

REM UUtHtt*ttttt*****t*tttHtt**t**ttt**t*i***t***tt******** 


G CALCULATION 


PHI 2 
FW 



2350 RETURN 

2398 REM 

2399 REM 

2400 REM 

2401 REM 

2402 REM 
2410 T*= .25: REM 
2420 W=2 : REM 
2430 RETURN 

2498 REM 

2499 REM 

2500 REM MATERIAL DEFORMATION PROPERTIES 

2501 REM 

2502 REM 

2505 PRINT : PRINT** MATERIAL : 21-6-9 SS" 

'2510 E=2 . 04E+O7 :REM YOUNG'S MODULUS 

2520 BN=110011 ! 

2530 M=- . 07 6 
2540 PP= .0383 
2550 N=5.75 
2960 RETURN 

2998 REM 

2999 REM 

3000 REM J-dA CURVE INPUT 

3001 REM 

3002 REM 

3060 ALP A* 7.226487E-14 
3070 BETA*-6 . 2973 138D-10 
3080 GMMA= .0000040561366# 

3150 RETURN 

3198 REM 

3199 REM 

3200 REM CRACK GEOMETRY 

3201 REM 

3202 REM 

3220 PRINT: INPUT "Initial A/T" ;ADT 

3230 IF (ADT>1) OR (ADT*1) THEN GOTO 3220 . 

3240 PRINT: INPUT "Initial A/C"; ADC 
3270 A=ADT*T 
3280 C«A/ADC 
3290 RETURN 

3998 REM 

3999 REM 

4000 REM SET-UP THE SCREEN 

4001 REM 

4002 REM 

4010 PI=3. 14159265359# 

4020 PIDT=PI/2 

4030 FOR DUM = 0 TO NUM 

4040 X (DUM) *C*COS ( DUM*PIDT/NUM) 

4050 Y (DUM) *A*SIN (DUM*PIDT/NUM) 

4060 XO(DUM)«X(DUM) : YO ( DUM) «Y (DUM) 

4070 NEXT DUM 

4075 Y(0)=0! : X (NUM) =0 ! : YO ( 0) *0 : XO (NUM) *0 
4170 SCREEN 9 

4180 WINDOW (-1,-. 25)— (3.8, 1.25) 

4190 LINE (0 , - . 25) - (0, 1 . 25) 

4200 LINE (-l,0)-(0,0) : LINE (2 . 6*C/T, 0) - (3 . 8 , 0) 

4210 LINE (-1,1)-(3.8,1) 

4290 RETURN 


PLATE DIMENSIONS 


PLATE THICKNESS 
PLATE WIDTH 



5 CLS: CLEAR: KEY OFF: SCREEN 0 

10 DEFINT D 

11 DEFDBL A-C, E-Z 

12 DIM X ( 50 ) ,Y(50) , NX ( 50 ) , NY ( 50) , ADEL (50) ,TH ( 50) , NXT ( 50) , NYT ( 50 ) 

13 DIM G(50) ,J(50) ,C(50) , EXT ( 50) , START ( 50) , INCA (50) ,YD(50) 

14 DIM XO ( 50 ) , YO ( 50 ) ,XDS(50) , JT(50) , XT( 50) , YT ( 50) , JDA(50) , JDAM(50) 


15 REM 

35 GOSUB 2000: REM MAXIMUM VPL AND VPL INCREMENT 

40 GOSUB 2 3 00: REM NUMBER OF POINTS ALONG THE CRACK FRONT 

45 GOSUB 2400: REM PLATE DIMENSIONS 

50 GOSUB 2500: REM MATERIAL DEFORMATION PROPERTIES 

55 GOSUB 3000: REM J-R CURVE INPUT 

60 GOSUB 3 2 00: REM CRACK GEOMETRY 

65 GOSUB 4000: REM SCREEN SET-UP 

7 0 REM 

71 REM .MAIN LOOP 

7 2 REM 

75 GOSUB 6000: REM PLOT THE CRACK FRONT 


81 REM 

82 FOR DUM=0 TO NUM 

83 J (DUM) =0 

84 NEXT DUM 

85 GOSUB 14000: REM 

86 REM 

90 GOSUB 7000 : REM 

91 REM 

92 GOSUB 10000: REM 

93 REM 

94 GOSUB 12000:REM 

95 REM 

96 GOSUB 8000: REM 

97 REM 

98 GOSUB 8500: REM 

99 REM 

100 GOSUB 15000 : REM 


101 REM 

105 GOSUB 9000: REM 

106 REM 

110 GOSUB 9500 : REM 


INCREMENT THE PLASTIC DISPLACEMENT 
DIRECTION OF CRACK EXTENSION 
CALCULATION OF EFFECTIVE C ALONG FRONT 
CALCULATE ELLIPTICAL ANGLES ALONG THE FRONT 
CALCULATE G ALONG THE CRACK FRONT 
CALCULATE Jtot ALONG THE CRACK FRONT 
CALCULATE dJ/da ALONG FRONT 
DELTA-A 

ADVANCE THE CRACK FRONT 


111 REM 
150 GOTO 75 

500 END 

501 REM ************************************************************** 

1998 REM 

1999 REM 

2000 REM VPL AND VPL INCREMENT 

2001 REM 

2002 REM 

2160 PRINT : INPUT "Enter the Maximum Plastic Displacement (in) ";VPLM 
2170 IVPL=VPLM/50 
2283 RETURN 

2298 REM 

2299 REM 

2300 REM NUMBER OF POINTS ALONG CRACK FRONT 

2301 REM 

2302 REM 

2303 NUM=10 

2304 NUM=NUM- 1 : REM THIS ACCOUNTS FOR NUMBERING BEGINNING AT 0 

2305 REM 

2335 TH ( 0) =0 : TH (NUM) *3 . 14 1592 65359#/ 2 



Appendix to DW Boatwright’s Thesis 


APPENDIX 


COMPUTER PROGRAM LISTING 



APPENDIX 
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4.4 Preliminary Results ol Profiling 

Two of the profiles produced has been included as figures 
(4.1a and 4.1b) for the sake of discussion. Figure (4.1a) is a 
profile of a crack taken at the surface. The precrack and the 
monotonic fracture regions are marked on the figure. Figure 
(4.1b) shows a profile of the same crack taken at the center of 
the cross-section. Again, the precrack, and the monotonic 
fracture regions are shown. The character of the two profiles is 
quite different. The monotonic fracture region at the surface 
(figure 4.1a) is at an angle to the precrack region, while that 
angle is not obvious for the central section (figure 4.1b). To 
illustrate a second observation, the net displacement has been 




Position, Crack Growth Direction, i, Inches Position, Crock Growth Direction, z, inchea 


Figure 4.1: Typical Crack Face Separation Profiles; (a) at surface, (b) at central cross-section. 
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Figure 4.2: Typical Crack Separation Plots; (a) at surface, (b) at central cross-section. 

calculated as the difference between the y-values of the upper 
and the lower parts of the profile, and figures (4.2a) and (4.2b) 
have been included below to register this second observation. 

In figures (4.2a) and (4.2b), the precrack region and the 
monotonic fracture region are again marked. The two separation 
profiles look dramatically different: a substantial region of 

stretch exists at the end of the precrack region and the 
beginning of the monotonic fracture region at the surface. This 
stretch is not apparent in the central cross-section. 

At this time, the profile data is being developed and 
observations are being made. 
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4.5 Future Efforts 

Profiling continues, and analysis of the profiles and 
synthesis of an approach to account for the constraint in the 
crack growth of these "planar" specimens. 



4.10 


4.7 References 


[ 4 . 1 ] H. A. Ernst, "Material Resistance and Instability 

beyond J-Controlled Crack Growth", Elastic 
Plastic Fracture: 2nd Symposium, Volume I— 
Inelastic Crack Analysis, ASTM STP 803, Shih and 
Gudas, Ed., (1983), Pp. 1-191-213. 

[4 2] Ernst, H. A., "Further Developments on the 

Modified J-Integral", ASTM STP 995, (1989), Pp. 

306-19 . 

[4.3] Rice, J. R. , "A Path Independent Integral and the 
Approximate Analysis of Strain Concentration by 
Notches and Cracks", Journal of Applied Mechanics , 
Transactions of the ASME, June 1968. 

[4.4] Tada, H. , Paris, P. C., and Irwin, G. R. , The _ 
Stress Analysis of Cracks Handbook, Second Edition 
(1985), Del Research Corp., St. Louis, MO. 

[4.5] Hutchinson, J. W. , "Plastic Stress and Strain 
Fields at the Crack Tip”, Journal of the Mechanics 
and Physics of Solids, Vol. 16 (1968), Pp. 337-47. 

[4.6] Rice, J. R., and Rosengren, G. F., "Plane Strain 
Deformation Near a Crack Tip in a Power-Law 
Hardening Material", Journal of the Mechanics and 
Physics of Solids, Vol. 16 (1968), Pp. 1-12. 



